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Thermal  non-equilibrium  processes  in  partially  ionized  plasmas  can  be  most  accurately  modeled 
by  collisional-radiative  kinetics.  This  level  of  detail  is  required  for  an  accurate  prediction  of  the 
plasma.  However,  the  resultant  system  of  equations  can  be  prohibitively  large,  making  multi¬ 
dimensional  and  unsteady  simulations  of  non-equilibrium  radiating  plasma  particularly 
challenging.  In  this  paper,  we  present  a  scheme  for  model  reduction  of  the  collisional-radiative 
kinetics,  by  combining  energy  levels  into  groups  and  deriving  the  corresponding  macroscopic  rates 
for  all  transitions.  Although  level-grouping  is  a  standard  approach  to  this  type  of  problem,  we 
provide  here  a  mechanism  for  achieving  higher-order  accuracy  by  accounting  for  the  level 
distribution  within  a  group.  The  accuracy  and  benefits  of  the  scheme  are  demonstrated  for  the 
generic  case  of  atomic  hydrogen  by  comparison  with  the  complete  solution  of  the  master  rate 
equations  and  other  methods.  ©2013  AIP  Publishing  LLC.  [http://dx.doi.org/10.1063/L4849417] 


I.  INTRODUCTION 

The  ability  to  model  plasma  flows  with  non-equilibrium 
chemistry  plays  an  important  role  in  a  number  of  applica¬ 
tions  including  but  not  limited  to  plasma  propulsion,1 
high-speed  reentry  flows,-  plasma-assisted  combustion  and 
interpretation  of  laser  diagnostics.1  In  order  to  better  under¬ 
stand  the  physical  characteristics  of  the  flow  and  the  cou¬ 
pling  with  chemistry  under  the  conditions  of  interests,  one 
needs  to  accurately  model  all  the  non-equilibrium  processes 
associated  with  the  atoms  and  molecules  (e.g.,  excitation, 
ionization,  and  dissociation)  through  collisional  and  radiative 
interactions.4-6  The  most  accurate  treatment  for  these  non¬ 
equilibrium  plasmas  requires  a  state-to-state  approach,7-13 
also  referred  to  as  collisional-radiative  (CR)  models,  in 
which  deviations  from  the  equilibrium  distribution  of  the  in¬ 
ternal  states  can  be  captured. 

These  CR  models,  although  very  accurate  from  a 
physics  point  of  view,  can  be  computationally  very  expen¬ 
sive  due  to  the  large  number  of  internal  states  for  which  the 
number  densities  must  be  computed.  For  an  atomic  plasma, 
these  states  correspond  to  all  the  electronic  excitation  levels 
of  the  various  neutral  and  ion  species  considered.  For  a  mo¬ 
lecular  plasma,  additional  degrees  of  freedom  such  as  the 
rotational  and  vibrational  modes  further  increase  the  level  of 
complexity.  In  addition,  these  molecular  degrees  of  freedom 
are  strongly  coupled  to  the  chemical  reactions.  For  example, 
vibrational  excitation  facilitates  dissociation  or  other  endo¬ 
thermic  reactions,  and  recombination  can  also  favor  the 
production  of  excited  states.  These  models,  derived  from  ab 
initio  cross  section  databases  for  all  elementary  processes, 
can  be  applied  to  a  wide  range  of  plasma  conditions  and 
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offer  more  complete  insight  into  the  non-equilibrium  effects. 
For  example,  a  recent  study  of  ionizing  shocks  in  Argon  by 
Kapper  and  Cambier10'11  demonstrated  that  this  level  of 
detail  is  needed  for  an  accurate  prediction  of  high-speed 
flows.  In  addition,  the  unsteady  coupling  of  the  hydrodynam¬ 
ics  and  CR  kinetics  leads  to  physical  phenomena  which  can, 
in  turn,  provide  additional  information  useful  for  model 
validation  and/or  experimental  measurements  of  various 
parameters. 

Due  to  the  large  computational  workload  involved  in 
solving  the  CR  master  equations,  simulations  incorporating 
state-to-state  kinetics  have  only  been  limited  to  zero-  or  one- 
dimension  with  a  few  exceptions  of  two-dimensional  calcu¬ 
lations.1 1141:1  For  example,  the  run-time  for  solving  a  set  of 
rate  equations  for  the  CR  kinetics  of  atomic  hydrogen  scales 
as  the  cubic  power  of  the  size  of  the  atomic  state  distribution 
function  (ASDF)  when  an  implicit,  backward-Euler  method 
is  employed.  While  better  scaling  laws  could  be  obtained 
with  iterative  and  more  approximate  schemes,  their  accuracy 
and  stability  for  extremely  stiff  problems  are  still  an  issue. 
The  development  of  very  efficient  and  accurate  schemes  for 
CR  kinetics  is  still  an  ongoing  research  topic  which  will  be 
presented  elsewhere;  here,  we  discuss  a  different  approach, 
consisting  of  lowering  the  complexity  of  the  calculations  by 
developing  a  reduced-order  kinetic  model  suitable  for  multi¬ 
dimensional  flow  calculations  while  maintaining  a  sufficient 
level  of  detail  required  to  accurately  model  the  plasma. 
Several  mechanism  reduction  schemes  have  been  proposed 
in  the  literature  with  applications  to  various  types  of  kinetics. 
Colonna  et  a/.16  utilize  a  two-level  distribution  model  to 
study  nitrogen  dissociation  rates  in  recombining  flows,  in 
which  all  the  vibrational  levels  except  for  the  last  level  are 
modeled  by  a  single  energy  equation  with  an  assumption  of  a 
Boltzmann  distribution,  and  the  last  vibrational  level  is  mod¬ 
eled  using  state-to-state  kinetics  to  take  in  account  the  non¬ 
equilibrium  effects  of  the  upper  states.  Magin  et  al.17  have 
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developed  a  rovibrational  collisional  (RVC)  coarse-grain 
model  to  characterize  the  internal  energy  excitation  and  dis¬ 
sociation  processes  of  nitrogen  flow  behind  a  strong  shock 
wave.  The  coarse-grain  model  is  derived  by  lumping  the 
rovibrational  energy  levels  into  groups,  in  which  the 
population  is  described  by  a  uniform  distribution.  Guy 
et  al. 18  proposed  a  multi-internal -temperatures  models  for  a 
vibrationally  non-equilibrium  flow,  in  which  the  vibrational 
distribution  is  divided  into  two  or  three  groups,  each  with  its 
own  vibrational  temperature.  Liu  et  al.,19  on  the  other  hand, 
proposed  a  mechanism  reduction  to  CR  models  based  on  the 
multi-group  maximum  entropy  principle  with  the  constraints 
being  the  macroscopic  parameters. 

In  this  paper,  we  examine  several  different  level  group¬ 
ing  schemes  for  the  state-to-state  kinetics  of  atomic  elec¬ 
tronic  states.  The  first  approach  is  similar  to  that  of  Magin 
et  al}1  for  the  rovibrational  collisional  coarse-grain  model 
and  therefore  is  based  on  uniform  (U)  binning  of  the  levels. 
The  second  approach  here  consists  of  grouping  levels  into 
groups  with  an  assumed  Boltzmann  (B)  distribution,  allow¬ 
ing  a  higher-order  description  of  the  ASDF.  In  this  case,  the 
effective  excitation  temperatures  are  evolved  in  time  by  con¬ 
serving  a  set  of  moments  of  the  distribution  function;  the 
most  obvious  solution  is  to  solve  for  number  density  and 
energy,  similar  to  the  approach  by  Guys  et  a/.18  However, 
we  will  show  that  a  different  set  of  moment  variables  of  the 
same  order  should  be  used,  due  to  the  specific  nature  of  the 
ASDF. 

The  method  developed  here  can  be  applied  to  a  wide 
range  of  state-to-state  kinetics  models  including  the 
RVC i s  i 7  ancj  vibrational  '  collisional  (VC)  models  or  the 
electronic  collisional-radiative  model.  ’  ’  In  the  interest 

of  simplicity,  we  consider  here  the  CR  model  of  atomic 
hydrogen,  using  classical  models  for  the  level  energies  and 
rates;  the  actual  values  of  these  parameters  are  unimportant 
here,  as  long  as  the  structure  of  the  ASDF  is  representative 
of  the  actual  species,  notably  the  geometric  progression  of 
the  level  energies  of  the  ASDF  and  the  stiffness  ratio.  The 
level  grouping  techniques  are  applied  to  reduce  the  cost  of 
solving  the  full  master  equations  and  the  results  are  com¬ 
pared  with  the  reference  solution  computed  from  the  full 
master  equations. 

The  rest  of  the  paper  is  organized  as  follows:  we 
describe  the  state-to-state  kinetics  and  numerical  solution  in 
Sec.  II,  while  in  Sec.  Ill  we  describe  the  various  mechanism 
reduction  methods.  We  compare  the  results  of  the  reduced- 
order  models  with  the  full  set  of  master  equations  in  Sec.  IV 
and  examine  the  issue  of  energy  conservation  in  Sec.  V. 
Finally,  a  summary  is  given  in  Sec.  VI,  while  the  derivation 
of  the  kinetic  rates  used  in  this  study  is  given  in  Appendix. 

II.  COLLISIONAL-RADIATIVE  MODEL 
A.  Definitions  and  rates 

As  mentioned  above,  we  consider  here  the  ASDF  of 
atomic  hydrogen  coupled  to  electron  impact  excitation  and 
ionization,  and  the  reverse  processes  (respectively,  deexcita¬ 
tion  and  recombination),  as  well  as  the  radiative  rates  for 
line  transitions  in  an  optically  thin  approximation.  Radiative 


recombination  is  neglected  and  all  radiation  absorption  is 
ignored,  as  is  free-free  (Bremsstrahlung)  emission,  since  this 
does  not  directly  affect  the  atomic  level  populations.21  The 
atomic  states  of  the  hydrogen  atom  are  listed  as  a  function  of 
their  principal  quantum  number  (n)  only,  following  the  Bohr 
atomic  model;  the  splitting  of  states  with  respect  to  orbital 
and  spin  numbers  is  ignored,  and  all  states  have  a  degeneracy 
gn=2n2.  The  states  number  from  n  =  I  to  oc  and  we  con¬ 
sider  a  finite  number  of  states  n=  1,  ...,  M<  oo  before 
reaching  the  ionization  limit.22  In  this  simplified  model,  the 
energy  of  each  state  is  given  as  En  =  I h\  1  —  l/«2),  as  meas¬ 
ured  from  the  ground  state  (Et  =  0),  and  we  will  denote  by 
I n  =  / n(\/n2  —  1  / M 2)  ~  Ih/ii2  the  energy  required  for  ioni¬ 
zation  of  level  n. 

The  population  density  N„  is  the  number  of  atoms  per 
unit  volume  of  a  state  n.  For  a  single  bound-bound  transition 
between  states  n  and  m  (m  >  n)  induced  by  electron-impact 
collisions,  the  rate  of  change  of  the  population  density  is  of 
the  form 


-jf=-^Un)N^e  +  PUm)NmNe.  (1) 

Hereafter,  we  will  use  the  convention  of  indexing  the  rates 
with  the  final  state  on  the  left,  and  the  initial  state  on  the 
right,  i.e.,  (f\i).  The  first  term  on  the  right  of  Eq.  (1) 
describes  the  loss  due  to  excitation  from  level  n  to  m,  as  a 
result  of  collisions  between  free  electrons  (of  number  density 
Ne );  the  second  term  describes  the  gain  due  to  collisional 
deexcitation  from  the  state  m,  with  number  density  Nm.  Note 
that  for  the  same  transition  between  the  levels  n  and  m,  we 
also  have 


dN, 


-^=+«Un)NnNe-P{n\m)NmNe. 


(2) 


If  there  were  only  two  states  to  consider,  Eq.  (1)  would  be 
the  entire  rate  of  change  for  level  n,  but  since  all  transitions 
involving  the  state  n  must  be  counted,  the  rate  of  change  for 
excitation  and  deexcitation  alone  involves  summing  up  the 
right  hand  side  over  all  levels  m^n.  At  equilibrium 
(Boltzmann),  the  ratio  of  population  densities  is 


K, 

N*n 


=  B„m{Te)  = 


8m  e-AE„m/kTe 
gn 


(3) 


where  A Enm  =  Em  —  En  is  the  difference  in  level  energies. 
For  electron-impact  processes,  the  rates  a  and  />  in  Eqs.  (1) 
and  (2)  are  functions  of  Te  and  are  given  by  Eqs.  (A9a)  and 
(A9b)  in  Appendix.  For  low  values  of  the  energy  gaps 
between  levels  (A Enm/kT  <C  1),  both  forward  (a)  and  back¬ 
ward  {p)  rates  become  very  large.  This  leads  to  a  wide  range 
of  time  scales  as  the  number  of  levels  is  increased,  and  to  a 
considerable  stiffness  in  the  system  of  equations.  For  exam¬ 
ple,  Figure  1  demonstrates  the  increase  in  both  the  maximum 
eigenvalue  (inverse  time  scale)  and  the  spread  of  values,  i.e., 
stiffness,  as  the  plasma  evolves  as  function  of  time. 
Additionally,  Eqs.  (A9a)  and  (A9b)  show  that  the  system  is 
strongly  diagonally  dominant,  in  the  sense  that  transitions 
with  small  changes  in  quantum  number  (m  —  n  ~  1)  have  a 
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FIG.  1.  Spectrum  of  eigenvalues  of  the  CR  system  versus  time,  during  con¬ 
stant-'/).  plasma  evolution  from  a  low-temperature  ASDF  and  low  electron 
number  density;  as  excitation  and  ionization  proceed,  the  upper  states  and 
Ne  increase,  yielding  a  rapid  growth  of  the  characteristic  frequencies. 


higher  rate  than  those  with  m  —  n  1.  In  fact,  a  fairly  good 
approximation  here  could  be  to  consider  a  ladder  process, 
i.e.,  transitions  between  neighboring  states  only,  but  this  is 
not  necessarily  applicable  to  other  atomic  configurations, 
and  this  approximation  is  not  used  here. 

For  ionization  and  recombination  processes,  the  rate  of 
change  of  the  population  density  for  level  n  is 

-*\Mn)NnNe  +  Fh+)N+N].  (4) 


The  first  term  on  the  right  side  is  the  loss  due  to  ionization  of 
that  level  by  electron  collisions  ( Ne ),  while  the  second  term 
is  due  to  the  capture  by  an  ion  (N+)  of  a  free  electron  (one 
factor  of  Ne),  in  the  presence  of  a  second  electron  (leading  to 
an  N2  dependence),  required  for  energy  conservation.  The 
equilibrium  for  ionization  and  recombination  (Saha)  involves 
a  different  relation 


=  Sn(Te) 

gn 


\  ll2 


3/2 


p-In/kTe 


(5) 


where  g+  is  the  degeneracy  of  the  ion  ground  state  (for 
atomic  hydrogen,  g+=  1).  Thus,  we  cannot  assume  that  the 
equilibrium  values  are  the  same  for  both  excitation/deexci¬ 
tation  and  ionization/recombination  processes.  Usually,  we 
can  have  Boltzmann  equilibrium  (3)  without  Saha  equilib¬ 
rium,  but  hardly  the  reverse,  mostly  because  it  takes 
more  energy  to  ionize  than  to  excite;  for  the  upper  states 
close  to  the  ionization  limit  (n  ^>1),  the  difference  is  less 
significant. 

Only  the  radiative  transitions  between  atomic  levels 
(“line,”  or  “bound-bound”  emission)  are  considered  here. 
For  each  bound-bound  transition  m  ■  n  (m  >  ri),  we  have 


dNn 

dt 


mi 


dNm 

dt 


A(n\m)N m- 


(6) 

(7) 


B.  Master  equations 

Once  all  the  macroscopic  rates  are  obtained,  we  can  con¬ 
struct  the  master  equations  describing  the  collisional- 
radiative  kinetics  of  all  levels.  In  this  study,  we  consider 
atomic  hydrogen,  which  has  only  one  ion  state,  and  only 
electron  collisions,  which  allows  us  to  remove  the  super¬ 
script  e  in  the  rate  definition  hereafter.  The  rate  of  change  of 
the  population  density  of  a  level  n  is  thus  written  as 

^  =  ^  ]  M(m\n)NeNn  +  y  '  ft  {n\m)N eN m  +  y  ^  ^{n\m)Nm 

m>n  m>n  m>n 

^  ^  M(n\m)NeNm  ^  ^  P{m\n)^e^n  ~  ^ 

m<n  m<n  m<n 

~a(+\ n)NeN„  +  P(n\+^N+N2e.  (8) 

Similarly,  we  can  write  another  equation  for  the  rate  of 
change  of  the  population  density  of  the  ions  according  to  the 
rate  of  ionization  or  recombination 

^  =  E  «(+l n)NeNn  -  E  P(n\+)N+N2e.  (9) 

n  n 

Finally,  the  electron  density  is  related  to  the  ion  density  by 
the  charge  neutrality  condition 

Ne  =  EZ?AV  (10) 

1 

We  will  compute  the  time  evolution  of  a  uniform  plasma;  if 
we  assume  a  constant  temperature  bath,  the  conservation 
equations  above  constitute  a  complete  set,  but  for  constant- 
volume  conditions — with  time  variation  of  the  tempera¬ 
ture — there  is  also  conservation  equation  for  the  electron 
energy,  which  will  be  examined  in  more  detail  in  Sec.  V. 
The  task  of  deriving  a  reduced  model  for  the  CR  kinetics 
aims  at  modeling  the  shape  of  the  ASDF  at  a  lower  computa¬ 
tional  cost  compared  to  that  required  to  solve  the  full  master 
equations,  while  maintaining  sufficient  accuracy  to  capture 
the  non-equilibrium  effects.  The  most  natural  way  to  accom¬ 
plish  this  is  to  partition  the  excited  states  into  groups  or 
“bins,”  therefore  reducing  the  number  of  variables  in  the  sys¬ 
tem.  Various  assumptions  can  be  made  about  the  internal 
structure  of  each  group,  i.e.,  the  distribution  of  states  within 
the  groups,  and  various  approaches  to  solving  the  group- 
based  variables  can  be  devised. 

C.  Numerical  solution 

Examination  of  Eqs.  (8)  and  (9)  reveals  that  the  full  sys¬ 
tem  of  ODEs  can  be  written  in  the  following  form: 

^=-J„-Xp  +  E;K PlXq  {/>,<?}€  {l,...,M,l+}, 

<7 

(11) 

where  Xp  is  the  /r-element  of  the  vector  of  conserved 
variables — for  the  set  of  master  equations  (8),  Xp  =  Np — and 
JP,  Kw  are  matrices  built  from  summation  over  all  possible 
transitions  between  levels,  and  are  themselves  functions  of 
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Xp.  Thus,  the  source  term  can  generally  be  decomposed  into 
a  linear  and  non-linear  terms.  The  system  (11)  can  be  solved 
using  a  variety  of  techniques.  In  this  study,  we  have  used  a 
backward  Euler  scheme  to  avoid  the  stiffness  of  the  CR 
kinetics.  Expanding  the  general  system  (11) 

dXp  =  -(J,  +  dip)  ■  (Xp  +  dXp) 

+  ]T(Kw  +  dKpq)(Xq  +  dXq)  (12) 

<7 

and  retaining  lst-order  terms  only,  we  obtain 

Apl.  ■  dXr  ~  -J p  ■ Xp  +  Y,  '  X,  (13) 

<7 

with  the  Jacobian 

V  =  (1  +  J P)Spr  -  K pr  +  Xp  -J2Xq  (  OXr)  ' 

(14) 

For  this  implicit  method,  there  is  no  stability  restriction  on 
the  time  step.  For  consistency,  all  the  simulations  shown  in 
this  paper  utilized  a  constant  time  step.  The  same  solution 
methodology  is  applied  to  the  various  cases  of  level  group¬ 
ings,  where  now  some  of  the  conserved  variables  in  the  set 
[Xp]  are  summations  over  the  levels  within  the  groups/bins, 
while  in  the  general  case  of  non-isothermal  plasma,  it  also 
includes  the  electron  energy  Ee. 


III.  LEVEL  GROUPING  STRATEGIES 
A.  Uniform  grouping 

Consider  a  group  of  M  individual  levels 
i  =  {no,  abbreviated  as  i  £  n  and  denote  the 

group,  or  “bin”  number  by  n;  hereafter,  n,  m,.. .  are  the  group 
indices  and  i,  j,...  are  level  indices.  This  first  approach  to 
model  reduction  is  essentially  a  zeroth-order  approximation 
of  the  internal2,  distribution  function,  where  only  one 
moment  variable,  either  the  total  number  density  of  the 
group  or  the  total  excitation  energy  of  the  group,  is  required. 
The  traditional  choice  is  to  conserve  the  total  number  density 
of  the  group,  i.e.,  Afn  =  /V,.  Using  Eq.  (3),  a  Boltzmann 

approximation  of  the  internal  partition  function  Z„  is 
obtained  by24 


6«0  ,'g„ 


-A E,/T. 


<  N 

ien  JV"0 


(15) 


where  A E,  =  Ei  —  E„0  is  the  difference  in  energy  between 
the  level  i  and  the  first  level  of  the  group,  n0.  The  approxima¬ 
tion  of  a  group  with  uniform  internal  distribution  is  equiva¬ 
lent  to  having  a  characteristic  group  temperature  Tn 
approaching  infinity,  compared  to  the  total  energy  width  of 
the  group,  i.e., 


2 


n 


(16) 


where  g„  is  the  overall  group  degeneracy.  The  simplest 
model  therefore  consists  of  assuming  all  levels  within  the 
group  to  be  distributed  uniformly,  i.e.,  weighted  by  the  level 
degeneracy 


Ni  =  -Afn. 

8n 


(17) 


The  rate  equation  for  a  group  n  is  obtained  by  summing  the 
master  rate  equations  (8)  and  (9)  for  all  the  levels  i  within 
the  group,  and  utilizing  relation  (15) 


dAf„ 

dt 


=  -NeA fn 
+  NeAfm 
-Af„ 
-NeAf, 


EEfE*oio+EEfE%i 

m>n  i£n  m<n  i£n  °nj£m 


8j  , 


m<n  i£n  j£m 

gi 


m>n  i£n  j£m  ‘ 


EEfE^io 

m<n  i£n  &nj£m 


\  '  gi 

2^~a(+ io 


;sk  8" 


+  Af, 

+  AavJE% 


EEEf%) 

m>n  i£n  j£m 


i+) 


(18) 


Similarly  for  the  ion  state,  one  obtains 


dN 4 


u7r=^E-A/’«  E -^EfE^ 


dt 


8n 


i+) 


(19) 


The  terms  within  brackets  in  Eqs.  (18)  and  (19)  contain 
effective  rates  for  the  groups,  which  can  be  pre-computed. 
For  example,  in  the  first  term  on  the  right-hand-side  of 
Eq.  (18) 


l(m\n) 


EfE- 

i£n  * n  j£m 


'010 


is  an  effective  excitation  rate  from  group  n  to  group  m.  Note 
that  since  this  model  does  not  require  computing  an  excita¬ 
tion  temperature  Tn,  all  the  effective  transition  rates  between 
the  groups  can  be  expressed  as  a  function  of  the  kinetic  tem¬ 
perature  Te  only.  It  is  important  to  emphasize  that  the  group¬ 
ing  of  levels  is  applied  on  the  high  energy  states  only;  thus  in 
any  simulation  we  must  choose  a  number  of  low-energy, 
“resolved”  levels,  as  well  as  a  variable  number  of  groups 
combining  the  upper  states.  The  number  of  discrete  states, 
the  number  of  groups  and  their  widths  are  variable  parame¬ 
ters  of  the  model,  whether  we  use  uniform  binning  as  above, 
or  Boltzmann  internal  distributions,  discussed  below.  In 
order  to  bound  this  parameter  space  (optimization  is  beyond 
the  scope  of  the  present  work),  we  need  to  provide  a  refer¬ 
ence  solution,  such  that  the  population  density  of  each  level 
can  be  compared  to  the  one  reconstructed  from  the  assumed 
internal  distribution  within  each  group.  Figure  2  shows  the 
evolution  of  the  electron  density  computed  from  the  master 
equations.  This  test  corresponds  to  a  strong  ionization  regime 


Le,  Karagozian,  and  Cambier 


Phys.  Plasmas  20,  123304  (2013) 


123304-5 


FIG.  2.  Time  evolution  of  the  electron  number  density  using  different  total 
number  of  atomic  levels.  The  electron  temperature  is  set  at  3.0  eV. 


ASDF  need  to  be  conserved.  The  selection  of  these  variables, 
however,  can  be  arbitrary.  Guy  et  a/.18  conserved  the  total 
number  density  of  the  group  and  the  average  excitation 
energy;  these,  respectively,  correspond  to  zeroth-  and 
first-order  moment  variables,  and  would  appear  to  be  a  natu¬ 
ral  choice.  Consider  the  total  number  of  states  J\fn — defined 
in  Eq.  (15) — and  the  total  energy  within  the  bin 
£n  =  ]T  0,  EiNj,  for  which  we  can  write  conservation  equa¬ 
tions,  derived  from  Eq.  (8) 


dN„ 

dt 


~Ne  Mn 


AEi/T„ 


EE^ — E 

m>n  ien  n  j£m 


0  p—AEj/Tn 

+EE^E — E%) 


m<n  ien 


jem 


(20a) 


and  the  time  evolution  of  the  ASDF  shows  an  increasing 
population  of  the  higher  atomic  levels  while  the  electron 
density  grows  exponentially.  It  also  demonstrates  the  effect 
of  the  number  of  levels  included  in  the  simulation,  i.e.,  using 
a  fewer  number  of  atomic  states  has  an  impact  on  delaying 
the  onset  of  the  electron  avalanche.  This  indicates  that  ioni¬ 
zation  from  the  high-energy  states  is  an  important  process, 
and  therefore  the  evolution  of  the  upper  states  must  be  accu¬ 
rately  captured.  We  could  always  increase  the  size  of  the 
ASDF  to  obtain  higher  accuracy,  but  with  diminishing 
return;  ultimately,  the  time-resolution  of  interest  and  the  ac¬ 
curacy  threshold  dictate  the  number  of  levels  required  in  a 
simulation.  The  mapping  between  the  practical  requirements 
and  ASDF  size  is  not  a  straightforward  matter,  but  is  an  issue 
beyond  the  scope  of  this  work.  Convergence  studies 
with  respect  to  the  size  of  the  system  showed  that  beyond 
20  levels,  there  were  no  discernible  differences  in  the 
results — see  Figure  2.  Thus,  we  chose  our  reference  solution 
to  be  the  one  obtained  for  20  levels,  and  all  level-grouping 
models  investigated  here  will  be  based  on  this  extent  of  the 
ASDF. 

B.  Boltzmann  grouping — Number  and  energy 


d£„ 

~dt 


~Ne  Mn 


EE^V— E£«aoi 


AEj/Tfj 


jem 


0.„-bE,/T„ 

-EE^E — E£‘%) 


jem 


(20b) 


For  sake  of  brevity,  we  did  not  write  the  entire  list  of 
contributions  in  Eq.  (20),  which  can  be  easily  derived  from 
Eq.  (18)  by  generalizing  the  weighting  factors  gj/gn  to 
gje~^EilT" / Zn,  and  similarly  for  other  groups.  By  solving  for 
total  number  and  total  energy  of  each  group,  according  to 
Eqs.  (20a)  and  (20b),  we  can  guarantee  direct  conservation 
of  both  mass  (total  number  of  levels  A f„)  and  energy  (£„). 
However,  this  approach  presents  some  problems  in  determin¬ 
ing  the  internal  Boltzmann  temperature,  as  will  now  be 
shown.  First,  let  us  define  a  total  group/bin  energy  measured 
from  the  lower  bound,  i.e.,  A£„  =  'Yhienfii  —  E^Ng,  the 
rate  of  change  of  this  shifted  energy  is  still  given  by  the 
right-hand-side  of  Eq.  (20b).  We  can  then  write 


A£n  =  gj  A £,  =  AT„(AE)n,  (21) 

Sn0  ien 


Several  assumptions  can  be  made  regarding  a 
Boltzmann-like  structure  within  the  group.  Panesi  et  al. 12 
and  Magin  et  al.25  rely  on  the  assumption  that  the  population 
within  a  group  follows  a  Boltzmann  distribution  at  the  ki¬ 
netic  temperature,  i.e.,  in  this  case,  Tn  =  Te.  This  approach  is 
only  valid  if  the  rates  of  exchange  between  the  levels  within 
the  group  are  much  faster  than  the  exchange  rates  with  levels 
outside  the  group;  otherwise,  one  could  then  assume  that  the 
entire  ASDF  is  governed  by  Te  and  is  always  in  Boltzmann 
equilibrium.  The  validity  of  this  assumption  is  highly  ques¬ 
tionable  for  atomic  state  populations.26  Furthermore,  when 
different  collision  partners  must  be  considered,  the  kinetic 
temperature  can  be  either  that  of  the  heavy  particles  or  the 
electrons  (e.g.,  electron-impact  excitation  and  heavy  impact 
quenching);  in  this  case,  choosing  either  one  of  the  kinetic 
temperature  can  impact  on  the  results. 

In  order  to  accurately  describe  the  population  of  a  group 
with  a  Boltzmann  distribution,  two  moment  variables  of  the 


where 


(A E)n  =  E  Y,  gi  A Ei  e-^1-  =  T2  E  in(^)  (22) 


2 -is! 


is  the  average  group  energy  measured  from  the  first  internal 
level.  Similarly,  a  specific  heat  at  constant-volume  can  be 
determined,  i.e., 


c^)=^(AE)n  =  T- 


E &w) 


VA Ei/T. 


=  t: 


(a  E2)n  (A E)l 


-ml 

(23) 


Since  N „  and  £n  are  conserved  variables,  we  obtain  new 
values  at  each  time  level  (k)  and  in  order  to  compute  the 
Boltzmann  temperature  T,„  we  need  to  iterate  the  equation 
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(A E)n(rn)  +  cv(rn)srn 


(24) 


where  T*  is  the  running  iterated  value,  until  convergence 
(ST*  ~  0).  However,  the  slope  of  the  curve  (AE)n(Tn)  is 
extremely  flat  at  low  temperature,  i.e.,  Cv  — >  0.  In  fact, 
when  T„  — >  0,  to  the  leading  order  we  have:  J\fn  ~  N„0 
o(l  +  e),  (A E)n  ~  o(e)  and  Cv(Tn)  ~  o(e),  where  e  =  e~*£'lTn 
is  a  small  parameter.  Therefore,  during  the  iterations  ST* 
=  o(e) / o{()  and  arbitrary  temperature  solutions  can  be 
obtained.  Our  studies  showed  that  indeed,  numerical  instabil¬ 
ities  prevent  us  from  obtaining  satisfactory  solutions  in  many 
test  cases.  While  it  is  possible  to  introduce  limiters  to  prevent 
unphysical  or  improbable  values  and  stop  the  iteration  coun¬ 
ters,  this  is  not  a  satisfactory  solution  to  the  problem.  We 
should  also  emphasize  that  the  problem  occurs  when  Tn  is 
small,  which  does  not  imply  that  electronic  levels  are  unpopu¬ 
lated,  since  we  may  very  well  have  small  internal  group  tem¬ 
peratures  as  a  result  of  initial  conditions  or  running  iterations, 
but  non-negligible  overall  electronic  excitation  (J\f„  /  0).~7 


the  group,  which  allows  us  to  separate  the  variables,  one  of 
o(l)  and  the  other  of  o(e).  Clearly,  we  have  now 


5n0  ienr 


-A Ei/Tn 


using  Nt=^gie-^-. 


Z'„ 


(27) 


In  order  to  evaluate  the  new  temperature  from  the  two  con¬ 
served  variables,  we  iterate  on  ST*  using  a  form  similar  to 
Eq.  (25) 


K(K)  + 


dZ'n 

dT„ 


a r'[k) 

*t:  = 


N 


(*) 


(28) 


However,  it  is  easy  to  see  that  since  jjjZ'  =  ^Z,  we  obtain 
a  similar  equation  to  Eq.  (26): 


r*  2 


ST*  ~ 


Z'n(T:){AE)n(T*) 


M"{k) 

gno-Z'n{T*n ) 


N 


C.  Boltzmann  grouping — Partitioning 

In  the  approach  above,  we  are  dealing  with  two  reduced 
variables  Af„  and  £„  (or  A£„)  which  are  both  summations 
over  the  internal  levels.  An  alternative  may  consist  of  keep¬ 
ing  one  of  the  level  populations  as  a  variable.  Therefore,  we 
could  instead  choose  for  each  group  n  to  conserve  the  popu¬ 
lation  of  the  lowest  level  in  that  group  N„0  and  Af„,  whose 
evolution  is  given  by  a  form  similar  to  Eq.  (20a).  To  evaluate 
the  Boltzmann  temperature  of  the  group,  we  now  have  at 
time  step  ( k ),  from  Eq.  (15): 


np 

§n0 


J2s'e 

i€n 


■AEi/T„ 


N W 

n0 

&nQ 


Zn{T {nk)) 


so  that  in  order  to  evaluate  the  new  bin  temperature  we 
need  to  solve 


In  the  same  limit  T„  — >  0,  both  numerators  and  denominators 
are  of  o(e)  and  the  temperature  iterations  are  again  unstable; 
this  was  verified  through  extensive  tests  under  a  variety  of 
conditions  and  configurations.  To  avoid  this  systematic  nu¬ 
merical  problem,  we  must  consider  another  way  to  evaluate 
the  Boltzmann  temperature  inside  each  group. 

Consider  instead  the  following  expansion  of  the  parti¬ 
tion  function  near  the  mean  relative  energy  value 
A E„  =  ZJ2isn8i^Ei-  Defining  Si  =  AEj  —  AEn  as  the 
shifted  energy  gap,  we  have 


8  i/T„ 


=  e 


—AEn/Tn 


i  - 


S, 

T„  2  Tl 


'gne-^IT'll+odS^/n) 


(29) 


zm  +  (c^)sT:=^gn  0 


dTn 


N 


until  convergence.  Using  Eq.  (22),  this  leads  to 

T*2 


ST !  ~ 


Zn(T*)(AE)n(T*) 


U{k) 

4  >  /7n 


(25) 


(26) 


where,  again,  the  dependencies  on  temperature  have  been  ex¬ 
plicitly  written.  At  low  T,„  the  denominator  is  o(e)(  1  +  e) 
and  the  numerator  is  a  difference  between  two  terms  of 
o(l  +  e).  Therefore,  the  iterative  procedure  is  again  numeri¬ 
cally  unstable. 

To  attempt  to  alleviate  this  problem,  we  have  examined 
yet  another  approach:  for  each  group  n  we  conserve  the  pop¬ 
ulation  of  the  lowest  level  in  that  group  N„a  and  A f'n,  the  total 
population  of  the  remaining  upper  states  n'  of  that  group, 
such  that  n  =  no  U  This  is  an  effective  partitioning  within 


where  g„  is  the  total  degeneracy — see  Eq.  (16).  Therefore, 
up  to  second-order  in  the  approximate  ratio  of  the  bin  width 
to  the  temperature,  the  partition  function  can  be  approxi¬ 
mated  by  a  single  exponential  function  and  the  relation  (29) 
can  be  inverted.  If  we  use  the  (N„0 ,  J\f)  pair  of  conserved  var¬ 
iables,  we  have 


Af® 


N 


go 


=  Zn{TW)  ^  gne-^ . 


(30) 


However,  the  left-hand-side  of  Eq.  (30)  is  o(  1  +  e),  and  the 
right-hand-side  should  be  as  well.  To  see  that  this  is  the  case, 
consider  the  first  terms  in  the  expansion  of  Eq.  (29)28 


Zn(T„) 


e-AE„ITn 


goe-(AEo-AE.)/TM  +  gie-m-*E.)/TH  +  .  .  . 


Since  AE  ~  AE\  and  AEf]  =  0,  the  right-hand-side  is 
o(e)[o(l/e)  +  1 +  •••]— o(l  +  e).  Again,  this  is  not  a 
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desirable  situation,  since  the  evaluation  of  the  group  temper¬ 
ature  Tn  is  of  the  form  1  /ln(  1  +  e),  and  is  subject  to  signifi¬ 
cant  errors.  Furthermore,  by  computing  the  average  gap  A E 
from  the  lower-bound  of  the  energy  bin,  the  requirement 
(5)  <C  T„  may  be  hard  to  justify  at  low  group  temperature. 

Instead,  we  can  take  advantage  of  the  self-similar  struc¬ 
ture  of  the  atomic  spectrum  (exact  for  hydrogen,  approxi¬ 
mate  for  other  atoms)  and  the  fact  that  the  energy  gaps 
become  narrower  as  the  level  index  increases.  Thus,  let  us 
define  the  average  energy  counting  from  the  first  level  above 
the  lowest  level,  as  obtained  from  Z'n,  defined  in  Eq.  (27) 


risking  fatal  numerical  errors;  this  is  possible  only  by  sepa¬ 
rating  the  lowest  and  upper  levels  within  the  group,  i.e.,  by 
performing  a  sub-scale,  internal  partitioning  of  the  group.29 
This  is  the  approach  used  here  for  the  last  Boltzmann  (here¬ 
after  denoted  as  B5)  group  we  investigated,  for  which  the 
appropriate  pair  of  conserved  variables  to  use  is  therefore 
(. N„0,Af'n ).  Note  that  it  is  also  possible  to  improve  on  the 
temperature  evaluation  by  incorporating  all  higher-order 
terms  into  the  definition  of  the  total  degeneracy,  i.e., 

Z'M  =  g'n(T„)e~^T" 


Z[,  =  ^gie-AE-^  =  e-M-P-Ygie-W-. 


iGn' 
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(34) 


By  definition  of  the  mean,  the  first-order  term  in  the  expan¬ 
sion  of  the  exponential  on  the  right-hand-side  should  be: 
Yhi&n'  8i&'i  —  0,  where  now  =  A E,  —  A E'n.  This  yields 

AE'n=~jJ2S'AEi  with  Sn  =  J2S'-  (32) 

iGn'  i>no 


Therefore,  A E’  differs  from  A E  only  by  a  normalization  fac¬ 
tor,  since  AEp  =  0.  Note  that  A E'  >  AEi  and  to  lowest- order, 
Z'(Tn)  ~  gfne~&E' ^T"  —  o(e).  Using  the  conserved  pair 
(. N„0,J\T '),  the  group  temperature  is  now  estimated  by 


A 


N, 


8n0  —  Z„(T„) 


jW  ^ _ 


A  E' 


1 


In 


N  n  gfl  o 
g'n  Nno 


ln(e) ' 
(33) 


This  is  now  a  stable  computation  when  e  — >  0.  Furthermore, 
the  approximation  (6)  <C  Tn  is  more  justifiable  since  the 
largest  value  (So  =  E„0  —  A E)  is  removed  from  the  average. 

We  see  that  we  now  have  the  means  to  compute  the  in¬ 
ternal  group  temperature  from  conserved  variables  without 


If  T*  is  the  running  iteration,  first  evaluated  by  Eq.  (33),  suc¬ 
cessive  estimates  of  T^'1  are  obtained,  using  Eq.  (34),  from: 


yW  _  t* 
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where 


(35) 


This  iterative  procedure  can  rapidly  converge  (as  demon¬ 
strated  in  our  tests)  because  we  have  an  excellent  approxima¬ 
tion  of  the  initial  temperature  from  the  lowest-order  direct 
evaluation  (33),  and  the  o(e)  term  has  been  factored  as  the 
leading  term  in  the  expansion.  In  other  words,  g'„(T„ )  is  a 
smooth  function  of  temperature  with  a  non-vanishing  gradi¬ 
ent,  allowing  gradient-descent  iterations. 

D.  Boltzmann  grouping — Effective  rates 

As  before,  the  master  equations  are  used  to  derive  the 
conservation  equations  for  the  two  new  variables  (N„0,Afn), 
by  setting  i  =  no  for  the  first  one,  and  summing  over  all  levels 
j  £  n!  in  the  second  case.  The  latter  yields  the  following: 
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(36) 


Note  that  we  have  used  the  total  number  Afm  =  Nmo  +  Af'm 
and  the  group  total  partition  function  Zm  =  gmo  +  Z'm  in  the 
expressions  on  the  right  hand  side,  only  as  a  way  to  group 
terms  and  lead  to  simpler  expressions;  the  conserved 


variables  remain  Nmo  and  Af'm.  Equation  (36)  takes  in 
account  all  the  interactions  between  the  groups,  assuming 
the  Boltzmann  distribution  approximation  within  each  group. 
The  effective  rates  for  group  transitions  can  be  expressed 
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(and  tabulated)  as  a  function  of  two  temperatures:  the  kinetic 
temperature  Te  and  the  group  excitation  temperature  T„. 
Notice  also  that  because  of  the  bin-averaging,  the  effective 
radiative  transition  rates  have  also  become  temperature- 
dependent  ( Tn ). 

Similarly,  the  rate  of  change  of  the  number  density  of 
the  ground  state  of  each  group  is 


dK, 

dt 


-=-NeN„  o 


Ne  J\fn 


EE  EE  Al«o) 

m>n  jem  m<n  jem 


+EE 


EE- 

m<n  j£m 

gje-^n. 


m>n  jem 


a(«ol /') 


-N„, 


E  Ea-m 

m<n  j£m 


+  A f 


EE 


8je 


m>n  j£m 


-A Ej/Tm 
"5  A(no\i) 


+NeK 


0,-AEf/Tn  0,„-AEi/Tn 

ElV-/?wo  +  E?V^wo 

ien1  n  ien1  « 


-AW,0[a(+|„o)]  +N;N+  [/?(„„!+)]  ■ 


(37) 


Again,  using  the  total  number  of  levels  J\fm  =  N„,n  +  J\f'm  on 
the  right-hand-side  allows  us  to  consider  together  transitions 
between  lowest  states  at  the  boundaries  of  different  groups 
(N„0  —Nm),  as  well  as  the  transitions  with  the  excited  sub¬ 
partitions  (7V„0  —  Afm)  and  simply  the  expressions.  Since  the 
ion  is  conserved  here  as  an  individual  state,  the  rate  of 
change  of  its  number  density  remains  the  same  but  can  be 
rewritten  in  terms  of  the  group  number  densities 
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(38) 


Each  term  in  brackets  in  Eqs.  (36)-(38)  is  an  effective  rate  for 
transfer  between  the  group  variables  (Nno,J\f'n),  Mn.  As  men¬ 
tioned  in  Sec.  Ill  A,  both  individual  levels  and  groups  (uni¬ 
form  or  Boltzmann)  are  considered  when  solving  the  ASDF. 
The  few  individual  states  are  the  lowest  in  the  energy  scale, 
with  the  largest  successive  gaps,  while  the  multitude  of  upper 
levels  is  distributed  into  a  variable  number  of  groups.  This  is 
justified  on  the  basis  of  the  kinetic  rates  (see  the  stiffness 
ratios  of  Figure  1),  and  as  justification  of  the  expansion  (29). 


TABLE  I.  Summary  of  level-grouping  models  investigated. 


Model 

Variables 

Equations 

Tn  evaluation 

U 

M„ 

(18)  and  (19) 

none 

B1 

(Mn,Sn) 

(20a),  (20b),  and  (38) 

Cv — unstable 

B2 

(N„0,M„) 

(37)  and  (20a) 

Cv — unstable 

B3 

(w»„X) 

(37)  and  (36) 

Cv — unstable 

B4 

(N„0,M„) 

(37)  and  (20a) 

Eq.  (29) — unstable 

B5 

(Nn0,M'„) 

(37)  and  (36) 

Eq.  (31) — stable 

able  to  provide  stable  and  satisfactory  solutions  for  all  test 
cases  was  model  B5,  using  a  sub-partition  of  the  group  into 
the  ground  level  n0  and  the  remainder,  and  the  use  of  the 
form  (31)  for  the  partition  function,  which  allowed  us  to  fac¬ 
torize  out  the  vanishingly  small  terms  at  low  T„.  Therefore, 
considerations  of  the  “equation  of  state”  of  the  Boltzmann 
group  dictated  the  correct  approach  to  use  here,  and  while  all 
the  models  explored  are  listed  in  Table  I,  only  the 
zeroth-order  uniform  binning  described  in  Sec.  Ill  A  and  the 
B5  models  are  shown  here  and  compared  to  the  reference  so¬ 
lution  obtained  from  solving  the  full  master  equations;  these 
are  indicated  as  (U)  and  (B)  models,  respectively. 

We  conducted  a  large  number  of  additional  tests  but  for 
the  sake  of  brevity,  we  are  showing  here  the  results  of  three 
representative  cases:  the  initial  conditions  are  summarized  in 
Table  II.  For  all  the  results  shown  in  this  section,  a  constant 
time  step  of  10-7  s  had  been  used  for  the  test  cases  in  the  ion¬ 
ization  regime  (cases  1  and  3),  and  a  time  step  of  10  5  s  was 
used  for  the  recombination  regime  (case  2);  the  same 
backward-Euler  scheme  of  Sec.  II C  was  used  throughout. 

As  indicated  in  Sec.  Ill,  the  reference  solution  is  based 
on  the  detailed  kinetics  for  20  atomic  levels,  while  the 
group-based  solutions  will  be  based  on  a  few  low  energy  lev¬ 
els  individually  monitored,  and  with  partitioning  of  the 
remaining  upper  states  into  a  variable  number  of  groups.  The 
first  test  case  is  the  iso-thermal  relaxation  in  the  excitation 
and  ionization  regime,  i.e.,  the  initial  population  of  excited 
states  and  electron  density  is  well  below  equilibrium.30  This 
test  case  is  the  same  as  the  one  shown  in  Figure  2  for  a  vari¬ 
able  number  of  electronic  levels,  solving  for  the  full  master 
equations  (8)  and  (9).  As  the  plasma  relaxes  towards  equilib¬ 
rium,  an  increasing  number  of  electronic  levels  become 
populated  and  the  electron  number  density  grows  exponen¬ 
tially,  until  an  ionization  cascade  occurs.  The  rates  increase 
very  rapidly  just  before  equilibrium,  and  the  system  becomes 
very  stiff,  as  shown  by  the  large  spread  of  eigenvalues  in 
Figure  1. 


IV.  ACCURACY  OF  UNIFORM  AND  BOLTZMANN 
METHODS 

A.  Isothermal  ionization  test  case 

In  the  previous  section,  we  have  discussed  several 
approaches  to  the  level  grouping  strategy;  these  are  summar¬ 
ized  in  Table  I.  This  sequence  of  models  was  developed  as  a 
result  of  preliminary  tests  and  the  failure  to  obtain  converged 
solutions  for  the  group  Boltzmann  temperature  T„  in  many 
instances.  Thus,  we  found  that  the  only  model  which  was 


TABLE  II.  Initial  conditions  of  test  cases.  For  all  cases,  the  total  atomic  den¬ 
sity  Nh  is  1021  m\ 


Case 

T 

1  e 

Xe  =  N+IN„ 

N„ 

i 

3  eV — isothermal 

10“9 

(1  ~xe)N„ 

to -2°nh 

for  n  =  1 

otherwise 

2 

1  eV — isothermal 

Saha  (3  eV) 

Boltzmann  (3  eV) 

3 

3  eV — isochoric 

10~9 

(1  -xe)NH 
to -20nh 

for  n  =  1 

otherwise 
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Figure  3  shows  a  comparison  of  the  number  densities  of 
all  the  atomic  states  for  the  iso-thermal  test  case  (1).  In  this 
simulation,  the  ground  state  and  the  first  4  excited  states 

(1, _ ,5)  are  conserved  as  discrete  levels  while  the  remaining 

upper  states  (6,..., 20)  are  partitioned  into  two  groups,  each 
of  which  has  either  a  uniform  or  Boltzmann  distribution. 
There  are  both  significant  and  subtle  differences  in  the  traces 


(a)  Full  solution  with  20  levels. 


(b)  Solution  with  5  levels  and  2  Uniform  groups 


time  (microsec) 


(c)  Solution  with  5  levels  and  2  Boltzmann  groups. 


FIG.  3.  Comparison  of  the  time  evolution  of  the  excited  states  during  the 
isothermal  heating  test  case  (Te  —  3  eV).  From  top  to  bottom:  (a)  full  solution 
with  20  levels;  (b)  solution  with  5  levels  and  2  uniform  groups;  (c)  solution 
with  5  levels  and  2  Boltzmann  groups.  The  first  excited  state — H(2) — is  the 
top  curve,  followed  by  the  next  higher  level,  etc. 


FIG.  4.  Comparison  between  the  solution  obtained  using  both  level  grouping 
approaches.  The  solid  line  represents  the  full  solution. 


of  the  upper  states.  First,  comparison  of  the  uniform  (Figure 
3(b))  and  Boltzmann  (Figure  3(c))  grouping  shows  the  influ¬ 
ence  of  the  assumed  internal  distribution,  as  the  recon¬ 
structed  levels  of  the  groups  are  clearly  separated  in  the 
uniform  case.  Second,  comparison  with  the  reference  solu¬ 
tion  of  Figure  3(a)  shows  that  the  Boltzmann  groups  are 
clearly  more  accurate.  Slight  differences  remain  in  the  very 
early  stages  of  evolution1'  below  1  /is  for  example  and 
around  10  /rs. 

The  combined  effect  of  the  number  of  resolved  lower 
levels  and  grouping  strategy  is  shown  in  Figure  4.  Generally 
speaking,  one  can  clearly  observe  a  dramatic  improvement, 
for  the  same  number  of  resolved  levels,  by  switching  from  a 
uniform  to  Boltzmann  group.32  By  selecting  the  time  of  max¬ 
imum  rate  of  growth  of  the  electron  density  as  the  approxi¬ 
mate  location  of  the  avalanche  ionization,  we  can  measure 
the  relative  error  in  density.  As  shown  in  Table  III,  the  error 
can  be  very  substantial  unless  there  is  sufficient  resolution  of 
the  ASDF  kinetics,  through  the  number  of  resolved  lower 
levels  and  a  higher-order  (B)  description  of  the  groups.  This 
is  important  when  comparing,  for  example,  with  time-gated 
experimental  results. 

By  conserving  more  discrete  states  and  reducing  the  size 
of  the  upper  state  groups,  the  results  are  of  course  signifi¬ 
cantly  improved.  This  is  to  be  expected  for  ASDF  kinetics, 
since  the  energy  gaps  are  larger  for  the  first  levels,  and 
grouping  together  these  states  would  be  less  accurate,  first  by 
yielding  excessive  bin  energy  widths  compared  to  mean 
energy  and  temperature  scale — violating  the  validity  condi¬ 
tion  for  the  expansion  (29) — and  also  by  disallowing 


TABLE  III.  Relative  error  on  electron  density  at  peak  rate  of  growth 
(approximately  33  /is). 


Method 

Error  (%) 

3  levels  +  1  U-group 

2618 

3  levels  +  1  B -group 

89.2 

4  levels  +  1  U-group 

165.8 

4  levels  +  1  B -group 

23.7 

6  levels  +  1  U-group 

20.9 

6  levels  +  1  B -group 

0.9 
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potential  deviations  from  Boltzmann  equilibrium  in  the  most 
populated  range  of  excited  states. 

There  can  of  course  be  variations  in  the  grouping  strat¬ 
egy,  but  in  all  cases  the  general  guidelines  of  keeping  the 
widths  of  the  groups  small  and  the  levels  with  the  largest 
energy  gaps  as  individual  states  are  perfectly  consistent  with 
the  objective  of  computational  cost  reduction,  since  the  dis¬ 
crete  lower  energy  states  evolve  more  slowly  and  the  upper 
states  are  numerous  and  have  similar  energy.” 

The  relative  accuracy  of  the  grouping  approaches  can 
also  be  seen  in  Figure  5  where  the  ASDF  is  plotted  at  four 
different  instances  of  time  corresponding  to  t=  10,  20,  30, 
and  40  ps.  Both  solutions  with  level  grouping  are  obtained 
from  using  3  atomic  levels  and  1  group  of  upper  states.  It  is 
clearly  seen  that  the  Boltzmann  group  gives  a  more  accurate 
representation  of  the  upper  states  distribution  during  the 
heating  process.  We  also  showed  in  Figure  5  the  results  of  a 
simplified  model  where  it  is  assumed  (see  Sec.  Ill  B)  that  all 
groups  have  the  same  internal  temperature,  equal  to  the  ki¬ 
netic  temperature,  i.e.,  T),(n)  =  Te.  "in  (dashed  line).  This 
assumption  is  clearly  violated,  as  shown  in  Figure  6, 
although  the  difference  remains  mostly  confined  to  the  upper 
states  distribution.  We  should  point  out  again  that  significant 
differences  would  be  expected  in  a  two-temperature  kinetic 
system,  i.e.,  including  heavy-particle  collisions. 

We  note  also  that  the  ASDF  from  the  full  solution  indi¬ 
cates  that  the  high  lying  states,  starting  from  the  third  excited 
state,  behave  like  a  continuum  state,  although  there  appears 
to  be  two  distinct  sub-groups  among  the  upper  states,  as  can 
be  seen  most  clearly  at  t  =  10  /ts.  This  suggests  that  the  upper 
states  are  most  effectively  resolved  by  two  groups  or  more, 
again  confirming  that  relatively  small  widths  of  the  groups 
are  preferable,  albeit  at  an  increased  computational  expense. 
Figure  6  further  illustrates  this  point  by  showing  the  evolu¬ 
tion  of  the  Boltzmann  temperatures  of  the  upper  states,  using 
here  4  discrete  atomic  states  and  partitioning  the  upper  states 
into  3  groups.  While  the  Boltzmann  temperatures  of  the  first 
two  groups  are  close  to  each  other,  the  temperature  of  the 
third  group  is  slightly  higher.  This  again  confirms  that  the 


FIG.  5.  The  internal  states  population  during  the  heating  process  at  various 
times.  The  solid  symbols  are  the  full  solution;  the  solid  lines  are  the  level 
grouping  with  Boltzmann  distribution;  the  doted  lines  are  for  level  grouping 
with  uniform  distribution;  dashed  lines  are  for  a  simplified  model  with 
T„  =  Te. 


0  5  10  15  20  25  30  35  40  45  50 

time  (microsec) 

FIG.  6.  Boltzmann  temperature  of  the  upper  states. 

upper  states  need  to  be  resolved  by  at  least  2  groups.  When 
the  system  is  near  equilibrium,  both  approaches  give  similar 
results. 

In  these  simulations,  we  have  assumed  that  the  plasma  is 
optically  thin  to  all  the  radiation  from  the  line  transitions. 
Spectral  signatures  being  a  major  diagnostic  tool  for  determin¬ 
ing  plasma  conditions,  it  is  important  to  know  the  CR  kinetics 
in  detail  in  order  to  match  experimental  data.  Usually,  this  is 
accomplished  by  post-processing  the  numerical  solution  with 
a  highly  resolved  spectral  code — including  radiation  transport 
(RT)  if  necessary — with  detailed  computation  of  line  shapes. 
This  approach  is  accurate  if  the  key  parameters  of  such  a  spec¬ 
tral  model,  in  particular  Ne  and  Te,  are  also  very  accurate.  As 
discussed  above  and  shown  in  Table  III,  our  Boltzmann 
grouping  procedure  provides  a  significant  improvement  over 
conventional  approaches,  leading  to  a  potentially  much  more 
accurate  spectral  signature  prediction  in  transient  and  non¬ 
equilibrium  plasma  conditions.  In  addition,  the  ASDF  solution 
is  much  closer  to  the  true  physical  state,  which  may  also  lead 
to  faster  integration  of  the  detailed  CR  kinetics  with  RT. 
These  will  be  investigated  in  the  future. 

Accurate  evaluation  of  the  radiative  emission  is  also  im¬ 
portant  during  the  computation  of  flow  dynamics,  from  simple 
reasons  of  power  coupling,  e.g.,  radiative  cooling.  Figure  7 
shows  the  radiative  losses  due  to  bound-bound  radiation 
from  the  upper  states  (5,..., 20)  to  the  first  three  atomic  states 
(1,  2,  3)  computed  by  grouping  all  the  upper  states  together  as 
a  single  group  with  a  Boltzmann  distribution.  Although  this  is 
a  somewhat  coarse  approximation  to  the  ASDF,  it  is  clear  that 
the  grouping  scheme  provides  an  excellent  approximation  to 
the  radiative  power.  An  accurate  reproduction  of  the  radiative 
spectrum  depends  inevitably  on  the  reconstructed  population 
of  the  atomic  levels  and,  as  can  be  seen  by  comparing  the  pro¬ 
files  in  Figure  3,  the  agreement  can  be  excellent. 

B.  Isothermal  recombination  test  case 

In  this  case,  we  performed  a  cooling  test  where  the 
plasma  is  suddenly  brought  down  from  3  eV  to  1  eV.  Thus, 
the  simulation  was  run  at  a  constant  temperature  ( Te  =  1  eV), 
while  the  initial  conditions  are  the  Boltzmann  and  Saha  equi¬ 
librium  values  at  3  eV;  these  are  exactly  the  conditions 
which  would  be  obtained  at  the  end  of  the  first  test  case  in 
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FIG.  7.  Radiative  loss  due  to  bound-bound  radiation  from  the  upper  states  to 
the  first  3  atomic  states.  The  lines  indicate  the  solution  obtained  from  the 
full  CR  kinetics.  The  dots  represent  solution  obtained  with  level  grouping 
(5  levels  +  1  group). 

the  absence  of  radiative  losses.  For  all  the  simulations  shown 
in  this  case,  a  constant  time  step  of  1 0  5  s  has  been  used. 

In  this  case,  the  cooling  process  occurs  very  rapidly  and 
the  plasma  is  in  a  deexcitation  and  recombination  regime; 
the  ground  state  and  the  electron  number  densities  are 
quickly  adjusted  to  their  new  equilibrium  values,  as  can  be 
seen  in  Figure  8.  Strictly  speaking,  since  bound-bound  radia¬ 
tion  is  assumed  to  be  optically  thin,  the  system  cannot  reach 
equilibrium.  However,  a  quasi-equilibrium  state  is  achieved 
at  approximately  1  ms,  after  which  the  bound-bound  radia¬ 
tion  is  the  dominant  net  rate  of  change  and  the  system  con¬ 
tinues  to  cool  down  at  the  radiative  time  scales.  Note  also 
that  the  uniform  grouping  is  significantly  less  time-accurate 
than  the  Boltzmann  method,  as  was  already  the  case  in  the 
ionization  regime — see  Figure  4. 

Figure  9  shows  the  evolution  of  the  excited  states  as 
function  of  time  for  reference,  uniform  groups  and 
Boltzmann  groups.  Once  again,  there  is  a  noticeable  discrep¬ 
ancy  between  the  reference  solution  and  the  uniform  bin 
model,  especially  concerning  the  red  curve  which  crosses 
other  levels  during  the  relaxation  process.  This  curve  is  the 
density  of  H(2),  the  first  excited  state,  and  is  an  effect  of  the 


FIG.  8.  Comparison  of  the  time  evolution  of  the  ground  state  and  the  free 
electrons  during  the  isothermal  cooling  process  (3  eV  — >  1  eV)  using  level 
grouping  with  uniform  and  Boltzmann  distribution  (3  levels  +  1  group). 


(a)  Full  solution  with  20  levels. 


(b)  Solution  with  3  levels  and  1  Uniform  group 


(c)  Solution  with  3  levels  and  1  Boltzmann  group. 

FIG.  9.  Comparison  of  the  time  evolution  of  the  excited  states  during  the 
isothennal  cooling  test  case  ( Te  =  1  eV).  From  top  to  bottom:  (a)  full 
solution  with  20  levels;  (b)  solution  with  3  levels  and  1  uniform  group; 
(c)  solution  with  3  levels  and  1  Boltzmann  group.  H(3) — is  the  bottom 
curve,  followed  by  the  next  higher  level,  etc.;  the  non-conforming  red  curve 
is  H(2).  The  higher  population  densities  as  the  level  index  increases  include 
the  increase  of  the  level  degeneracy. 

strong  radiative  decay  of  this  state.  Notice  that  the  plot  starts 
at  t=  10“ 5  s,  i.e.,  the  first  implicit  time  step,  but  already  the 
solution  is  far  from  the  Boltzmann  equilibrium  which  is  the 
initial  condition  at  t  =  0,  such  that  there  is  a  population  inver¬ 
sion  with  respect  to  H(2)  for  many  upper  states.  Notice  also 
that  the  time  scale  is  logarithmic,  and  the  processes  consider¬ 
ably  slow  down  as  the  electron  density  drops  significantly. 
Because  we  are  considering  only  electron  impact  collisions, 
the  ASDF  essentially  becomes  “frozen”  in  a  quasi-static  but 
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FIG.  10.  Comparison  of  the  time  evolution  of  the  excited  states  number 
densities  during  the  isothermal  cooling  process.  The  lines  are  full  solution. 
The  dots  represent  solution  obtained  with  level  grouping  (3  levels  +  1  B 
group). 


non-equilibrium  state.  If  collisions  by  heavy  particles  were 
also  considered,  these  would  rapidly  become  the  dominant 
process,  leading  to  faster  relaxation  towards  equilibrium. 
However,  in  some  case  of  rapid  plasma  expansion,  similar 
“frozen-in”  non-equilibrium  distribution  functions  of  the 
ASDF  could  be  obtained. 

To  better  appreciate  the  accuracy  of  the  Boltzmann 
grouping  procedure.  Figure  10  shows  the  evolution  of  sev¬ 
eral  excited  states  compared  to  the  exact  solution  and  simi¬ 
larly  to  the  “heating”  (ionization)  case,  excellent  agreement 
was  obtained.  In  this  simulation,  the  first  3  atomic  states 
(1,  2,  3)  are  conserved  as  discrete  levels  and  the  upper  states 
(4,. .  .,20)  are  lumped  into  1  Boltzmann  group. 

Finally,  we  show  in  Figure  11  the  snapshots  of  the 
ASDF  during  the  recombination.  Contrary  to  the  case  of  ion¬ 
ization,  the  upper  states  are  not  depleted  but  enhanced 
instead — as  expected,  since  the  recombination  proceeds  pref¬ 
erentially  onto  the  upper  states.  As  a  reflection  of  the  obser¬ 
vation  made  for  Figure  10,  the  agreement  is  excellent  for  all 
atomic  states. 


V.  ENERGY  CONSERVATION 
A.  Effective  rates 

The  systems  of  equations  ( 18) — (19)  and  (36)-(38) 
describe  the  complete  evolution  of  the  ASDF  but  for  an  iso¬ 
thermal  plasma.  In  the  more  general  case,  the  ASDF  kinetics 
are  coupled  to  the  energy  of  the  system;  here,  this  includes 
only  the  total  energy  of  the  free  electrons  Ee.  Thus  for 
constant-volume  or  constant-pressure  conditions,  there  must 
be  an  evolution  equation  for  the  energy  or  enthalpy  (only 
constant-volume  kinetics  are  considered  here).  We  must  then 
exert  care  that  the  formulation  exactly  conserves  energy,  i.e., 
that  +  P2n  £pp  at  any  time  level  ( k )  remains  the  same 
within  numerical  round-off  errors.  If  we  were  dealing  with 
only  electron-impact  collisions,  it  would  be  sufficient  to  sum 
the  energies  of  all  levels  using  the  new  population  densities 
at  the  end  of  the  time  step,  compute  the  difference  and  assign 
the  change  to  Ee.  However,  there  are  two  obvious  problems 
with  this  scenario:  (a)  when  other  collision  partners  must  be 
accounted  for,  or  when  the  electrons  themselves  are  parti¬ 
tioned  (e.g.,  for  non-Maxwellian  kinetics),  one  must  be  able 
to  correctly  apportion  the  changes  in  energy,  e.g.,  to  Ee  and 
Eh  (for  heavy  particles)  and  (b)  for  large  time  steps,  there  is 
no  guarantee  that  the  subsequent  change  in  Ee  is  physically 
acceptable,  i.e.,  Ep'1  =  Ep~ 1 :  +  5Ee  >  0.  We  must  therefore 
include  an  evolution  equation  for  Ee  (and  another  for  Eh  if 
heavy  particle  collisions  are  included),  which  must  then  be 
fully  coupled,  so  that  the  Jacobian  of  the  system  (14) 
includes  derivatives  of  the  rates  with  respect  to  Ee,  through 
the  variation  of  Te. 

Energy  conservation  can  be  satisfied  if  the  construction 
of  the  source  term  on  the  right-hand-side  of  the  master  equa¬ 
tions  also  satisfies  it.  Thus,  we  must  explicitly  construct  the 
energy  source  term  from  the  master  equations,  as  was  al¬ 
ready  described  briefly  in  Eq.  (20).  The  same  procedure  is 
used,  with  the  understanding  that 

dEe  v  -  d£n 

dt  2—/  c if 

n 


Thus  we  can  combine  contributions  as  follows: 


dEe 

—  Ne  J\f  n 

dt 


+X!  AEj‘P(j\i) 


J2  ^j^W) 

n 


-AEi/T„ 


j€m 


-AEj/T, 


z, 


j&n 


(39) 


FIG.  11.  Snapshots  of  the  ASDF  at  several  times  during  the  cooling  process. 
The  dots  represent  the  full  solution.  The  solid  lines  are  the  solution  obtained 
using  the  level  grouping  with  Boltzmann  distribution.  The  broken  lines  are 
the  solution  obtained  using  the  level  grouping  with  uniform  distribution. 


where  A E);  =  Ej  —  £).  Note  that  in  the  case  of  excitation 
from  level  |  i),  i.e.,  the  first  summation  in  Eq.  (39),  A  Ej,  >  0, 
while  A Eji  <  0  in  the  second  term  for  de-excitations  from 
that  level.  We  can  then  construct  another  set  of  effective 
rates,  this  time  for  the  energy  equation.  Using  the  sub¬ 
partitioning  of  model  B5,  the  rates  derived  from  the  first 
term  on  the  right  of  Eq.  (36)  are 
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~P(m'\n!) 


0.-AEi/T„ 

E?V-EAV<#> 


jCm' 


p.p—AEj/T,, 


yem' 


(40a) 


(40b) 


These  rates  enter  the  evolution  equation  for  Ee  as 

=  -NeKY,  “fm'|n')  ~  N‘’YY,  4  ■  (41) 


Note  that  the  same  formulation  applies  for  uniform  groups 
by  taking  the  limit  T„  — »  00,  and  summing  over  the  complete 
set  n={no,n'}.  The  rate  of  energy  change  can  also  be 
expressed  as 


^( m'\n ')  ^(m'\nr)  '  ^{m!\n!)  i  (4^) 

where  0E(m/u/)  is  of  course  given  by  the  effective  rate  for  the 
conserved  number  densities 


~  _  Sie  SE"^T"  _ 

a(m'|n')  —  Y  Y  Z'  ^ 

j£m!  i£nr  n 


Equation  (42)  defines  an  average  energy  „/),  transferred 
during  excitation  of  levels  of  group  n!  to  levels  of  group  m! , 
which  can  be  tabulated  as  function  of  the  initial  T„  and  colli- 
sional  ( Te )  temperatures.  This  approach  was  successfully 
used,  for  example,  for  vibrational  non-equilibrium.  ’4 

B.  Isochoric  ionization  test  case 

The  third  test  case  of  Table  II  was  designed  to  test  for 
energy  conservation.  In  this  case,  the  energy  loss  and  gain 
due  to  collisional  processes  are  taken  into  account  in  the 
conservation  equation  for  the  electron  energy.  The  evolu¬ 
tion  now  proceeds  at  constant  volume,  and  the  electron 
temperature  changes  rapidly,  as  seen  in  Figure  12.  The  ini¬ 
tial  conditions  are  the  same  as  those  of  the  first  test  case, 
and  the  system  is  initially  far  below  Boltzmann  and  Saha 
equilibrium  However,  contrary  to  the  isothermal  case,  the 


FIG.  12.  Ne,  Te  evolution  in  constant-volume  case. 
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FIG.  13.  Cumulative  and  instantaneous  relative  errors  in  energy  conserva¬ 
tion — test  case  3. 

initial  excitation  and  ionization  processes  deplete  the  elec¬ 
tron  energy  and  the  system  “freezes”  rapidly,  and  the 
excited  states  remain  at  a  low  population  density.  If  an 
external  heating  source  was  applied  (e.g..  Ohmic  heating), 
the  system  would  more  closely  resemble  the  isothermal  test 
case,  and  the  system  would  become  stiff  again.  Here,  we 
are  mostly  concerned  with  testing  energy  conservation  and 
to  simplify  the  analysis,  the  radiative  rates  were  removed 
from  the  kinetics,  so  that  no  radiative  energy  losses  were 
present. 

We  can  monitor  the  error  by  comparing  the  values  of  Ee 
at  the  end  of  each  time  step  with  the  total  potential  energy 
contained  in  the  electronic  states,  by  reconstruction  of  the 
level  populations.  Figure  13  shows  both  the  accumulated 
error  (symbols)  and  the  one  at  each  time  step  (solid  line); 
this  test  was  conducted  with  5  resolved  levels  and  2 
Boltzmann  groups,  and  all  computations  were  performed 
with  a  constant  time  step  of  10-10s,  using  the  same 
backward-Euler  integration  scheme.  ’5 

The  error  is  certainly  acceptable,  but  it  is  not  commen¬ 
surate  with  numerical  round-off,  which  we  would  have 
expected  if  the  scheme  was  exactly  energy-conserving.  By 
comparison,  the  cumulative  error  in  energy  was  below  Hr  11 
when  solving  the  full  master  equations  without  level 
grouping. 

While  the  exact  solution  consists  of  summing-up  the 
contributions  from  each  individual  level,  leading  to  the  rate 
of  change  expressed  by  Eq.  (20b).  However,  we  are  not  using 
here  the  internal  energy  £n  as  a  conserved  variable,  and  we 
must  be  careful  that  the  procedure  be  consistent  with  our 
definition,  or  reconstruction  of  the  internal  energy.  The 
corrected  procedure  is  described  next. 

C.  Corrected  energy  rates 

Consider  for  example  the  change  in  electron  energy  due 
to  excitations  and  de-excitations,  and  let  us  examine  first  the 
case  of  uniform  grouping. 

r]P  _  _  _  _  _  p 

=  X]  ^(m\n)-^nNe  +  Y,  Y/  P (n\mY n,N e ■  (43) 

m>n  n  m>n  n 
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There  are  two  formulations  of  the  effective  rates  of  energy 
transfer: 


d(AE)n, 

dt 


(AE2) 


(A E)l 


dT„  _ 
dt 


dTn_ 

'  dt 


(48) 


Formulation  1: 


From  Eq.  (31), 


~aE 

aH») 

^EEf^-^Hi  10. 

jam  ian  °n 

(44a) 

ft  (m\n) 

=  -«>%)■ 
jam  ian  °m 

(44b) 

Formulation  2: 

m\n ) 

=  (Em  -  En  ^2  ~a(j\i)j 

jem  ien 

(45a) 

h  m\n) 

=  ( Em  -  £,)EE/%. 

jam  ian  °m 

(45b) 

where  E„  =  j-E,  and  similarly  for  Em. 

Only  the  second  formulation  is  exactly  energy- 
conserving.  This  is  quite  clear  because  in  that  case,  the  term 
on  the  right  side  of  Eq.  (45)  is  the  product  of  the  change  in 
number  density  of  the  groups  (dAf„  / dt)  and  the  difference  in 
average  group  energy  ( E ).  Energy  conservation  follows  from 
the  definition  of  the  total  group  energy  £n  =  EnAf„.  Thus, 
the  model  assumptions  constrain  us  to  choose  the  appropri¬ 
ate  formulation  of  the  effective  rates  for  energy  change  that 
is  consistent  with  the  definition  of  group  energy. 

Let  us  now  examine  the  case  of  the  Boltzmann  grouping 
(B5),  using  the  pair  of  conserved  variables  (N„0,Af'n)',  the 
rates  of  energy  exchange  must  therefore  be  consistent  with 
the  electronic  energy  defined  from  these  two  variables,  and 
with  the  equation  of  state  used  to  describe  the  internal  parti¬ 
tion  (i.e.,  T „).  We  start  with  the  conservation  of  the  group 
energy 


dZ'n  gno 

'dK 

KdN  no 

=  K 

A  En  d\ng'n 

dT„ 

dt  N„a 

dt 

N„0  dt 

T2n  dT„ 

dt 

(49) 


Inserting  into  Eq.  (47), 


kti  d[AE)^ 
11  dt 


C  tT2 


(a E'„  +  T2^1ji 


dT„ 


dM’ 


dt 


KdN njL 

Nno  dt 


(50) 


We  can  now  combine  with  the  other  terms  of  Eq.  (47)  to 
obtain  an  expression  which  only  depends  on  the  rates  of 
change  of  the  conserved  variables  (Nno,Af'n).  Defining 


L’=- 


>  rr2 

/7 


AE'n  +  n 


dlng'n 

dTn 


and 


_  .  A f n: 

N~ 

iyn0 


(51) 


and  adding  the  contribution  from  the  ground  state  of  the 
group,  we  obtain 


d£„ 

dt 


[Eno  -  av]  ^  +  [Eno  +  (A E) n,  +  &,]  ^  .  (52) 


One  can  then  identify  the  rates  of  change  of  the  population 
density  with  the  effective  rates.  Considering  transitions 
between  groups  n  and  m  >  n,  and  using  a  similar  expression 
for  d£m/dt,  we  have 


d£„ 

dt 


d 

dt 


(NnoEno+N'n{E)n) 


=  E„ 


dN„0 
1  dt 


(E)n, 


dK  ,  d(E)n, 


dt 


dt 


(46) 


Note  that  the  averaging  ()n,  is  done  for  the  remaining  levels 
above  the  ground  level  nQ  of  that  group.  We  can  write  a  simi¬ 
lar  equation  for  the  total  energy  measured  from  the  ground 
state  of  that  group,  i.e., 


clA£n 

dt 


(47) 


The  first  term  in  Eq.  (47)  describes  the  change  in  group 
energy  from  the  global  change  in  population  of  the  group, 
i.e.,  (E) ndJ\[ n / dt .  The  last  term  describes  the  change  of  the 
internal  structure  of  the  group  as  a  result  of  the  collisional 
transitions,  since 


*(m0|n0)  —  [Emo  mm'  En 0  +  «„']  •  Ot(m0\n0) 

=  £(m0|«o)  -  ®(m0|no)i  (53a) 

*(m'|n0)  =  [Em0  E  (A E)m,  +  —  Eno  +  COn/]  •  S(m/| „0) 

—  ^(m'\no)  ’  Ct(m'\no):  (53b) 

®(m0|«')  =  [Emo  —  <J>m'  ~  E„0  —  (AE) n,  —  £„']  •  S(mQ|„q 

—  £( mo\n ')  '  ^(mo\nr)i  (53c) 

®(m'|;i')  =  [Emo  (A E) m,  +  —  E„0  —  (A E) n,  —  (,„/]  •  ®(m'|n') 

=  ’  y-(m'\n')  •  (53d) 


It  is  instructive  to  examine  the  limit  of  infinite  Boltzmann 
temperatures;  in  this  case, 

Zn<  ->  g'„,  C £n>,0)n'  — >  0  and  {AE)n,^AE'„ 

and  similarly  for  m! .  Equation  (52)  becomes 
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(54) 


where  we  have  also  used  the  fact  that  in  that  limit, 
N n0/ 8 n0  =  Nri/ gn,  and  used  the  definition  of  the  average 
group  energy — see  Eq.  (45).  Since  a  similar  equation  is 
found  for  d£m/dt,  the  combination  exactly  yields  Eq.  (45). 
Thus,  we  have  verified  that  by  taking  the  limit  TniTm  >  oo, 
we  recover  the  uniform  group  model. 

For  ionizations  and  recombinations,  a  similar  procedure 
can  be  found.  Considering  the  change  in  electron  energy  due 
to  ionization  and  recombination  from  and  to  the  group  n,  we 
have 


where  /,  is  the  ionization  potential  for  level  i  and  (I)n,  is  the 
group  ionization  potential  averaged  over  the  sub-partition  n' . 
Using  /,  =  IH  —  Ej  =  /„„  —  A Eh  it  is  easy  to  see  that 

(/)„,=/„„-(  AE)n,  and  ^  =  -Cvj.  (56) 


Equations  (49)  and  (50)  are  still  valid,  and  using  again  the 
definitions  (51),  we  obtain  the  final  form 


dEe\ 

~dt  J 


\/no  T  ^4?'] 
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Note  the  similarity  with  Eq.  (52).  The  effective  rates  are 
therefore 

*(+|«o)  =  [^"0  +  °V]  ■  ®(+|«o)  =  HMno)  '  “(+l«o)’  (58a) 

®(+|«')  =  [I*  ~  i^E) n'  ~  “=«']  ■  ®(+| n')  —  £(+| n')  '  ®(+|n')- 

(58b) 


Examination  of  equations  (53)  and  (58)  reveals  that  the  over¬ 
all  procedure  consists  of  replacing  the  energy  of  the  group’s 
ground  state  «o  and  sub-partition  n!  by  effective  energies  for 
the  energy  exchange 

E„o=E„0-a>n'  and  E,,'  =  E„0  +  (AE)n,  +  (59) 

Thus,  the  effective  rates  of  energy  transfer  become 


*  (mo|  no) 

(Emo 

Ena) 

'  a(m0|H0)i 

(60a) 

a(n?'|n0) 

{Em, 

-En0) 

'  ®(m'|n0)i 

(60b) 

a(mo\n’)  ~ 

{Em o 

~  En') 

'  ®(m0|n')> 

(60c) 

to-10  1(T9  1CT8  1(T7  10~6  10-6 

time  (sec) 


FIG.  14.  Cumulative  and  instantaneous  relative  errors  in  energy  conserva¬ 
tion — test  case  3 — with  revised  formulation. 

m'\n ')  i^m1  (60d) 

and  for  ionization: 

S(+M  =  ~  Etio)  '  ^(+K)’  (61a) 

af+| =  ( Ih  ~  En' )  •  «(+!„')•  (61b) 

The  use  of  effective  group  energies36  provides  a  straightfor¬ 
ward  approach,  and  the  effective  rates  of  energy  transfer  for 
all  transitions  (including  de-excitations,  recombination,  and 
radiative  transitions)  can  now  be  expressed  in  a  simple  form. 
Note  that  Eq.  (61)  is  similar  to  the  case  of  uniform  grouping 
(45)  and  since  we  have  already  demonstrated  that  we  can 
recover  the  uniform  grouping  case  in  the  limit  of  infinite 
temperatures,  we  have  achieved  here  a  fully  consistent 
model. 

We  are  now  left  with  the  task  of  verifying  energy  con¬ 
servation  with  this  revised  approach.  Using  the  same  test 
case  (3),  we  now  find  a  much  smaller  level  of  error,  as  can 
be  seen  from  Figure  14 — compare  with  Figure  13 — that  is 
the  characteristic  of  the  level  of  numerical  round-off.  Note 
that  the  cumulative  error  sums  the  absolute  values  of  the 
stepwise  error  (Lt  norm),  and  is  therefore  a  maximum  bound. 
Figure  15  shows  the  effect  of  bin  size  on  the  relative  error; 


FIG.  15.  Cumulative  relative  errors  in  energy  conservation  as  function  of 
group  sizes;  revised  formulation. 
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this  observation  is  similar  to  the  one  made  regarding  the  ac¬ 
curacy  of  the  ASDF — see  Figure  5,  i.e.,  smaller  group  widths 
are  preferred.  However,  it  is  clear  that  even  for  one  or  two 
bins,  the  error  on  energy  conservation  remains  very  small. 

VI.  CONCLUDING  REMARKS 

In  this  paper,  we  described  a  model  reduction  mecha¬ 
nism  for  the  collisional-radiative  kinetics  of  the  ASDF,  by 
grouping  electronic  states  into  groups  and  deriving  the  cor¬ 
responding  macroscopic  rates  to  take  in  account  all  the 
transitions.  While  level-grouping  is  a  commonly  used  and 
necessary  procedure  when  dealing  with  a  very  large  num¬ 
ber  of  atomic  levels,  as  in  high-temperature  plasma,37'38  the 
procedure  is  commonly  based  on  uniform  grouping  (i.e., 
simple  average  over  the  level  degeneracies).  A  higher-order 
description  of  the  internal  structure  of  the  groups  was 
developed  here  by  assuming  a  Boltzmann  distribution  of 
the  levels  within  the  group,  with  different  temperatures  for 
each  group.  This  approach  was  shown  here  to  have  superior 
features  to  the  uniform  grouping  procedure,  as  summarized 
below. 

Numerical  instabilities  in  the  limit  of  low  internal 
temperature  led  us  to  the  design  of  this  non-traditional 
approach.  Instead  of  conserving  the  total  energy  of  a 
group,  we  used  a  novel  approach  using  conserved  varia¬ 
bles  from  a  sub-partitioning  of  the  group  between  the 
lowest  level  of  that  group  (N„0)  and  the  remaining  states 
(A f'n).  Combined  with  an  appropriate  expansion  of  the 
sub-group  partition  function,  this  model  was  able  to  very 
rapidly  determine  the  group  temperatures  in  a  stable  and 
accurate  fashion  for  all  conditions,  solving  a  problem  that 
is  particularly  vexing  for  atomic  collisional-radiative 
kinetics,  due  to  the  structure  of  the  energy  levels  of 
atomic  plasma. 

The  two  grouping  schemes  (uniform  and  Boltzmann) 
were  implemented  on  an  electronic  collisional-radiative 
model  for  atomic  hydrogen  and  the  results  were  compared 
with  the  full  solution  of  the  master  equation  on  a  large  num¬ 
ber  of  numerical  tests,  of  which  a  representative  sub-set  was 
shown  here.  Both  schemes  showed  agreement  with  the  solu¬ 
tion  of  the  master  equations  and  have  excellent  convergence 
properties.  However,  substantial  accuracy  improvements 
were  obtained  at  minimal  computational  cost  with  the 
Boltzmann  grouping.  Detailed  tests  of  energy  conservation 
revealed  the  need  for  a  revised  approach  to  the  construction 
of  effective  rates  of  energy  transfer.  This  was  made  neces¬ 
sary  because  the  total  group  energy  was  not  part  of  the  set  of 
conserved  variables,  and  consistency  requirements  led  us  to 
a  new  formulation  involving  effective  level  energies  and  av¬ 
erage  group  energies  between  which  energy  exchange 
occurs.  We  derived  this  new  formulation  and  showed  that  it 
was  simple  to  implement,  and  consistent  with  the  uniform 
grouping  method. 

For  the  case  of  atomic  hydrogen  studied  here,  it  was 
sufficient  to  mix  a  few  low-lying  states  with  two 
Boltzmann  groups  for  all  the  high  energy  states  and  be 
able  to  capture  the  correct  ionization  kinetics  for  all  the 
states  and  the  radiative  spectrum.  This  significantly 


reduces  the  computational  cost  associated  with  solving 
the  kinetics  and  therefore  can  be  applied  in  multidimen¬ 
sional  and  time-resolved  flow  calculations,  with  accurate 
coupling  to  radiative  processes.  The  reduction  schemes 
presented  in  this  work  could  also  be  applied  to  other  set 
of  kinetics,  e.g.,  rovibrational  collisional  and  vibrational 
kinetics,  although  in  this  case  the  levels  are  distributed 
much  more  uniformly  on  the  energy  axis,  in  which  case 
reduction  schemes  and  group  temperature  determination 
are  more  straightforward.  Future  work  includes  a  straight¬ 
forward  extension  of  the  approach  to  non-hydrogenic  and 
multi-stage  ionization,  as  well  as  the  application  of  more 
accurate  and  more  efficient  time  integration  schemes, 
currently  under  development.  Further  optimization  of  the 
scheme,  such  as  dynamic  re-partitioning,  is  also  under 
exploration. 


ACKNOWLEDGMENTS 

We  wish  to  acknowledge  the  support  of  the  Air 
Force  Office  of  Scientific  Research  (AFOSR),  Grant  No. 
12RZ06COR  (PM:  Dr.  F.  Fahroo)  for  this  work. 


APPENDIX:  COMPUTATION  OF  KINETIC  RATES 

The  classical  form  of  the  cross-section  for  energy 
exchange  between  a  free  electron  and  the  atom,6  leading  to 
an  excitation  from  level  n  to  level  m  >  n  is 

=  (4 na20)  '  (3f™),  (Al) 

nm t 


where  a0  is  the  Bohr  radius,  e  is  the  energy  of  the  free 
electron,  A Enm  =  Em  —  En  is  the  energy  gap  between  n  and 
m  and  fnm  is  the  oscillator  strength 


f 


nm 


32  1  1  1 


3nV3n5  m 3 


(A2) 


The  free  electrons  are  assumed  to  follow  an  isotropic 
Maxwellian  distribution  fe(£e) 

fe(se)dee  =  2  ,  sl/2e^Ee/kTedse,  (A3) 

MkTe)3/2 

where  me  is  the  electron  mass,  ee  =  mev2J 2  and  Te  is  the 
temperature.  The  rate  of  excitation  is  obtained  by  averaging 
over  the  distribution  function 


a 


( m\n ) 


»oo 

<tenm(E)V<f(Ve)dve, 

E„m 


(A4) 


leading  to 


aH«)  =  (4J“o)tV 
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(3f  nm  )«A  nm  i 


(A5) 


where 
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Ve  = 


(  8kTe 
\nme 


^ nm  = - E  i(xnm),  and 

X/im 


E  i(*)  = 


y 


-dy. 


(A6) 


Here,  ve  is  the  mean  thermal  electron  velocity,  xnm 
=  A Enm/kTe  and  Et  is  the  exponential  integral.  The  reverse 
rate  can  be  found  from  detailed  balance 


Re  = 
P(nH  m2  e 


•  aH«) 


(A7) 


Pm 


4  n2h 3 
4  aon  n6 

n  m2kTe 


(A13b) 


The  rates  of  radiative  transitions  between  levels  can  also  be 
obtained  classically  for  the  hydrogen  atom.5  The  spontane¬ 
ous  emission  rates  from  an  upper  level  m  are 


A(h|/7j) 


87t2e2\  g„ 

*  I  '  i 

.  mec 3  J  gm 


1.6  x  1010 

m3n(m2  —  n2) 


(A14) 


The  expression  on  the  right  is  for  atomic  hydrogen  only. 


We  use  the  low  temperature  approximation6  (xnm  3>  1 ) 


(AS) 


in  which  case 


aH«)  — 


Ana^ 


32 
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The  factor  in  brackets  is  an  upper  bound,  which  is  reached 
for  the  upper  states  when  xnm  — >  0.  Another  scale  is  the  fac¬ 
tor  1  a! kTt,  in  xnm,  which  is  effectively  responsible  for  the 
stiffness.  If  that  factor  is  very  low  (high  temperatures),  all 
rates  are  of  the  same  order;  at  low  temperatures,  the  expo¬ 
nential  term  dominates  and  the  range  of  time  scales  is 
increased. 

The  cross-section  for  ionization  by  electron  impact  has  a 
form  similar  to  Eq.  (Al),  i.e., 


a;  =  (4tm§)  ^  In) .  (A10) 

This  leads  to  an  ionization  rate  coefficient6 

«(+!„)  =  {4nc$)ve(j^J  \l/(xn).  (All) 

The  final  state  (+|  is  an  ionized  state,  i.e.,  where  one  electron 
initially  bound  to  the  atom  has  reached  the  ionization  limit 
(«  =  oo )  and  is  part  of  a  free  continuum  of  states.  Using  the 
principle  of  detailed  balance,  the  reverse  (recombination) 
rate  is 


Pm  - 


A  a2h 3 
7i  m2kTe 
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Using  the  same  low  temperature  approximation  (A8),  we 
obtain6 
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